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Certain types of topological superconductors and superfluids are known to host protected Ma- 
jorana zero modes in cores of Abrikosov vortices. When such vortices are arranged in a dense 
periodic lattice one expects zero modes from neighboring vortices to hybridize and form dispersing 
bands. Understanding the structure of these bands is essential for the schemes that aim to employ 
the zero modes in quantum computation applications and in studies of their strongly interacting 
phases. We investigate here the band formation phenomenon in two concrete models, describing 
a two dimensional p x + ip y superconductor and a superconducting surface of a three-dimensional 
strong topological insulator (Fu-Kane model), using a combination of analytical and numerical tech¬ 
niques. We find that the physics of the Majorana bands is well described by tight binding models 
of Majorana fermions coupled to a static Z 2 gauge field with a non-trivial gauge flux through each 
plaquette, in accord with expectations based on very general arguments. In the case of the Fu-Kane 
model we also find that, irrespective of the lattice geometry, the Majorana band becomes completely 
flat at the so called neutrality point (chemical potential coincident with the Dirac point) where the 
model exhibits an extra chiral symmetry. In this limit the low energy physics will be dominated by 
four-fermion interaction terms which are permitted by symmetries and may arise from the Coulomb 
interaction between the constituent electron degrees of freedom. 


I. INTRODUCTION 

Topological superconductors attract our attention in 
part because they often host unpaired Majorana zero 
modes DKZI- These in turn exhibit a number of intrigu¬ 
ing physical properties, including a possibility to encode 
quantum information in a way that is robust to environ¬ 
mental decoherence [8] as well as to perform a limited set 
of quantum gates that are topologically protected mn- 
Thus far experimental evidence for Majorana zero modes 
(MZMs) exists in quasi one-dimensional systems includ¬ 
ing semiconductor quantum wires [T2lfT8l and wires com¬ 
posed of magnetic adatoms on a superconducting surface 
mm- 2D heterostructures made of strong topological 
insulators (STI) and conventional superconductors (SC) 
are beginning to also show promise [2TH29] . The latter 
systems are predicted to host unpaired MZMs in cores 
of Abrikosov vortices na and the quantum information 
stored in the zero mode subspace can be manipulated by 
performing adiabatic exchanges - “braiding” - of the indi¬ 
vidual vortices. How exactly one performs such braiding 
operations, initializes the system and reads out the re¬ 
sulting quantum state remains largely an open question 
but is one of considerable interest. 

A necessary first step towards the long term goal of 
storing and manipulating quantum information in the 
Hilbert space spanned by MZMs bound to vortex cores 
is to understand and characterize these systems from the 
point of view of their electronic structure. With this goal 
in mind we investigate here the electronic structure of 
topological superconductors in the presence of Abrikosov 
vortex lattices, paying particular attention to the fate of 
the MZMs associated with individual vortices. Specifi¬ 
cally, we study two different models of electrons in a 2D 
topological superconductor in the presence of a periodic 


vortex lattice. One describes a 2D spin polarized p x -Yip y 
superconductor and the other a SC surface of a 3D STI, 
also known as the Fu-Kane model m • Although not 
explicitly studied here, we expect our results to apply 
to other realizations of 2D topological superconductors, 
such as those predicted to occur in heterostructures com¬ 
bining spin-orbit coupled semiconductors, ferromagnetic 
insulators and ordinary superconductors EDGE]. Be¬ 
cause of the interplay between the orbital effects of the 
applied magnetic field B that is necessary to establish the 
vortex lattice and the spatially varying phase field 6{r) of 
the SC order parameter this turns out to be a problem of 
considerable subtlety and complexity. A variant of this 
problem in the p x + ip y superconductor and an 5-wave 
SC with Rashba spin-orbit coupling, but neglecting the 
applied magnetic field, has been studied in recent works 
[35H35] . It is however well known [361 E7] that the vor¬ 
tex lattice is thermodynamically unstable in the absence 
of B. (It is analogous to a system of charged particles 
without an appropriate neutralizing background.) 

To understand the electronic structure of a physical 
vortex lattice the magnetic field must be properly in¬ 
cluded. We do this here using a technique of the singular 
gauge transformation [38] first developed in the context 
of high-T c cuprate superconductors with d-wave symme¬ 
try and subsequently applied to both s- and p-wave su¬ 
perconductors 39 . The chief advantage of this technique 
is that it treats the magnetic field B and the SC phase 
field 0(r) on an equal footing. Indeed we find results for 
low-energy Majorana modes that differ in several impor¬ 
tant aspects from the results of Refs. [33H35] obtained 
while neglecting the B field. Most notably the flat Ma¬ 
jorana bands predicted in Ref. [33] and the strong band 
structure anisotropies seen in Ref. [35] are not present in 
the full solution of the problem. We conclude that mag- 
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netic field must be included in any calculation that aims 
to correctly capture the physics of Majorana zero modes 
in a realistic, thermodynamically stable vortex lattice. 
We note that semiclassical treatment of a p x + ip y super¬ 
conductor has recently been carried out (including the 
B field) j40| and showed results consistent with our fully 
quantum mechanical calculations. 

Our main results can be summarized as follows. De¬ 
noting a Majorana zero mode operator associated with 
a vortex core positioned at Rj by 7 j we find that in the 
presence of the vortex lattice the physics of these modes 
is well described by a tight binding model of the form 

^kin = ^ ' (-0 

iJ 

Here Uj = UjSij can be decoposed into a real symmet¬ 
ric matrix Uj representing the hopping strength while 
Sij = = ±i are Z 2 gauge factors. The imaginary 

unity present in the latter is dictated by the Majorana 
commutation relations 

{7i,7j} = 2<%, 7i=7i, (2) 

and the requirement that Hkin be hermitian. The sign 
ambiguity in Sij arises from the fact that one can perform 
a local Z 2 gauge transformation 7 j —>• — 7 j without affect¬ 
ing the zero mode commutation algebra ([ 2 |. A product 
of factors along a closed path, however, represents 
a Z 2 gauge flux that is gauge invariant, and therefore 
in principle observable. It is fixed by the microscopic 
Hamiltonian and can be thought of as analogous to the 
magnetic flux expressed through Peierls factors in lattice 
models of charged particles. Our main finding is that for 
the MZMs in the vortex lattice the Z 2 gauge flux through 
a general polygon formed by n vortices is given by the 
Grosfeld-Stern rule [41 

4*ij = 2 " ( n — 2), (3) 

polygon 

previously derived for MZMs in the Moore-Read frac¬ 
tional quantum Hall state [42 , whose effective theory is 
analogous to the spin polarized p x + ip y superconduc¬ 
tor. Eq. (|3| indicates a non-zero Z 2 gauge flux | and 
7 r through an elementary triangular and square plaque- 
tte, respectively, of the Majorana lattice. This in turn 
implies that the gapped phases of the Hamiltonian 0 
are typically topologically non-trivial with the occupied 
bands characterized by a non-zero Chern number. 

We note that the Majorana tight-binding model Eq. 
0 has been previously derived [13 E] for vortex lat¬ 
tices present in the Kitaev spin model on the honeycomb 
lattice [45]. It has been conjectured in these works that 
similar results should apply to other systems supporting 
localized Majorana mode arrays, but this conjecture has 
not yet been verified. Our work shows that the tight- 
binding model Eq. 0 with the Z 2 gauge structure § 
does apply to Abrikosov lattices in p-wave superconduc¬ 
tors and the SC surfaces of topological insulators. 


When the applied magnetic field is well below the 
upper critical field H c 2 we furthermore find that the 
hopping amplitudes Uj are significant only between the 
first and second nearest neighbors (nn). More gener¬ 
ally, preserve all vortex lattice symmetries and ex¬ 
hibit an exponential decay ^ e~ dij ^ with the distance 
dij = \Ri~ Rj \ and £ the SC coherence length, superim¬ 
posed on the RKKY-type oscillation with a period close 
to the Fermi momentum kp of the underlying normal 
metal. As already mentioned we find no sign of flat bands 
resulting from a subset of vanishing Uj predicted in Ref. 
|33l or anisotropies that break the underlying vortex lat¬ 
tice symmetry predicted in Ref. [35]. These effects ap¬ 
pear to be artifacts introduced by an approximation that 
neglects the magnetic field B. They may be present in 
small clusters of vortices if such can be stabilized in the 
absence of B but are not characteristic of a physical vor¬ 
tex lattice that retains its stability in the thermodynamic 
limit. 

The results described above pertain to both the p x -\-ip y 
superconductor and the Fu-Kane model. The latter 
shows an additional interesting feature when tuned to 
the neutrality point, reached when the chemical potential 
fi of the STI coincides with the Dirac point of the sur¬ 
face state. As noted previously [461148] . the model then 
exhibits an additional “chiral” symmetry which changes 
the topological classification of its zero modes from Z 2 to 
Z. Physically, this means that in the Fu-Kane model at 
neutrality MZM hybridization is prohibited and the Ma¬ 
jorana band must remain flat irrespective of the geometry 
of the vortex lattice. We confirm by explicit numerical 
calculation that this is indeed the case. We also find 
that when slightly detuned from neutrality MZMs form 
a weakly dispersive band with a narrow bandwidth pro¬ 
portional to the chemical potential p measured relative 
to the Dirac point. The Majorana flat band obtained by 
tuning a single parameter constitutes an interesting sys¬ 
tem because, just like in the fractional quantum Hall liq¬ 
uids USED], the kinetic energy of the particles becomes 
quenched and the nature of the ground state is deter¬ 
mined by interactions or disorder effects. If the sample 
is sufficiently clean so that disorder can be neglected and 
when interactions are present the system is inherently 
strongly correlated. Some consequences of these strong 
interactions in various ID and 2D vortex lattice geome¬ 
tries have been explored in recent studies [5TH56] . Effects 
of disorder on the Majorana tight-binding model Eq. 0 
have also been studied by several groups and interesting 
disorder-induced phases have been found |57H6l]. 


II. MAJORANA ZERO MODES IN VORTEX 
LATTICES 

Majorana zero modes associated with the individ¬ 
ual vortices in p x + ip y superconductor and the Fu- 
Kane model have been amply discussed in the literature 
[3 [331 S3 HI]. In this section we give a brief overview 
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of some key results and then focus on the effect of the 
applied magnetic field on the collective behavior of the 
zero modes in such lattices. Using approximate analyti¬ 
cal techniques we show how the low-energy Hamiltonian 
([]]) emerges in this setting and give expressions for the 
overlap integrals Uj and the Z 2 gauge factors Sij valid in 
a physical vortex lattice that includes the magnetic field. 


A. Spin polarized p x + ip y superconductor 

This is the simplest model of a 2D topological su¬ 
perconductor possibly relevant to Sr 2 Ru 04 [62], the A 
phase of superfluid 3 He [63] and the Moore-Read frac¬ 
tional quantum Hall state 0. The system is described 
by a second quantized Hamiltonian 

u = JMtH(r)V r , V r = , (4) 

where c}. is a spinless fermion creation operator. The 
Bogoliubov-de Gennes (BdG) Hamiltonian has the form 



where h = p 2 /2m — p is the kinetic energy operator and 
A = kp 1 {A(r), d x + id y } is the p x + ip y pairing operator 
with A(r) the SC gap function. It respects the particle- 
hole symmetry generated by S = r x K where r are Pauli 
matrices in the Nambu space and K denotes complex 
conjugation (S 2 = 1). In addition it is invariant under 
the global U(l) transformation H —)> e zr x He~ %T x when 
accompanied by a phase shift A —>► Ae~ 2lx 

We are interested in the solutions of the BdG equation 

H$(r) = E<S>(r) (6) 

in the presence of vortices in the SC order parameter 
A(r) = | A(r)\e t6 ^ with 6{r) is the SC phase. In the 
presence of singly quantized vortices located at spatial 
positions [Rj\ we may write 

K r ) = E ^fc( r ) = ar § ( r “ R k)- ( 7 ) 

k 

In addition A (r) vanishes at the center of each vor¬ 
tex and can be well approximated E3 as A (r) ca 
A o Ylj fanh (\r — Rj |/£). When the vortices are well sep¬ 
arated so that the smallest distance d £ then we may 
look for the low energy solutions of the BdG equation 0 
separately in the vicinity of each vortex. To this end we 
approximate the phase field near vortex j as 

0(0 ^ <Pj(r) + @j ( 8 ) 

where ©j = Vk(Rj) is the phase contributed by 

all other vortices in the system. Since by definition this 
contribution varies slowly near Rj it is permissible to 


approximate it by a constant. For future reference we 
also note that in view of Eq. 0 we can write 

©i = <KRj) (9) 

if we define 0(Rj) as being evaluated slightly to the right 
of the actual vortex position Rj , thus avoiding the sin¬ 
gularity at the vortex center. 

For a vortex j the zero mode BdG wavefunction can 
thus be written as 03 

/ e *(^+ 0 U 2 - 7r / 4 ) \ 

$ j( r i) = /( r i) y e -i( Vj +e,:/2-*/4)J , ( 10 ) 


where Vj=r — Rj and 


kp 


fir) = d ^Ji{k F r)ex.p 


[ |A(r , )|dr / 
vf Jo 


( 11 ) 


The quasiparticle operator 7 j = f d 2 r^j(rj)^ r associ¬ 
ated with the zero mode has the property 7 J = 7 j and is 
therefore Majorana. At this level of approximation each 
vortex contains a single Majorana mode. These modes 
have zero energy, obey canonical commutation relations 
0 . and are separated from the rest of the spectrum by 
a minigap A m — A q/Ef- Under adiabatic exchange 
vortices exhibit non-Abelian exchange statistics charac¬ 
teristic of the Ising anyons mm- 

To understand the electronic structure of the zero 
modes beyond the independent vortex approximation 
we must consider non-vanishing overlaps between their 
wavefunctions (10). To leading order the resulting low- 
energy Hamiltonian takes the form of Eq. ([l]) with 


iij = (Siim-). 


( 12 ) 


Refs, m HE] studied these overlap amplitudes between 
two vortices in various limits and found characteristic os¬ 
cillatory RKKY-type behavior with an exponential decay 
as can be expected on the basis of Eq. (11). Biswas [33] . 


in addition pointed out a specific dependence on phase 
angles Qj of the form 


. . ©; 

Uj oc sin v 


©, 


(13) 


This follows directly from the spinor structure displayed 


in Eq. (10) and has important consequences for the collec¬ 
tive behavior of MZMs in situations with many vortices, 
such as in the Abrikosov lattice. 


B. Fu-Kane model 


Fu-Kane model describes a SC surface of a 3D STI and 
is defined by the second quantized Hamiltonian 0 with 
a BdG Hamiltonian of the form 


H FK (r) 


h 

A* 


A 

-<jyh*cry 


(14) 
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where cr are Pauli matrices acting in the physical spin 
space, h = vp • cr — p and A = diag(A(r), A(r)). 
The Hamiltonian acts on a four-component spinor = 
(c^ r , c^ r , cj r , — cj r ) T in the combined spin and Nambu 
space. It respects the particle-hole symmetry generated 
by S = r y a y K (S 2 = 1 ) as well as the global U(l) sym¬ 
metry defined below Eq. ©• 

The Hamiltonian ( |l4| ) is known to support unpaired 
Majorana zero modes in singly quantized vortices and 
antivortices m- Their general properties have been ex¬ 
plored in Ref. [48]. Here we focus on the regime close to 
the neutrality point p = 0 where Hyk exhibits an extra 
chiral symmetry generated by U = r z a z . As a result the 
structure of MZMs becomes particularly simple, 


/ 

= fo(rj) 


e i(e j /2-n/4) 

0 

0 


\ 


(15) 


y_ e -i(© i /2-7r/4) J 

with f(r) = Aexp[—^ f (Q r | A(r')\dr']. The chiral sym¬ 
metry has an important consequence that the overlap 
amplitudes t^ between distinct MZMs exactly vanish at 
the neutrality point [46, 08]. When the symmetry is 
weakly broken by a small non-zero p then the overlap 
becomes m 


Uj = ipFij sin 



(16) 


where the integral is taken along the straight line be¬ 
tween Ri and Rj. In the remainder of this Section we 
shall justify this replacement in greater detail and we also 
evaluate the gauge invariant factors ujij in some specific 
situations of interest. 

In solving this problem we follow the classic procedure 
originally developed by Peierls [64] to include the mag¬ 
netic field in the tight binding model for electrons moving 
in the ionic lattice. It relies on a key assumption that the 
applied magnetic field is sufficiently weak so that the vec¬ 
tor potential A(r) can be replaced by a constant A(Rj) 
for the purposes of obtaining the individual zero mode 
bound state 4>j(r). Such a constant vector potential can 
then be removed from the kinetic energy term in Eqs. 
© and ( p~4| ) by the gauge transformation if x(r) is 
chosen such that 


VxiW = (19) 

In this gauge the zero mode is an eigenstate of the same 
Hamiltonian as in the absence of B except the SC phase 
is now given by 

%) = - 2 Xj{r). ( 20 ) 

k 

As before, near vortex j we can separate the slowly vary¬ 
ing part of the phase field and approximate it as 


with Fij = J d 2 rf 0 (r - Ri)fo(r - Rj). 


6{r) ~ipA r ) + ®i ( 21 ) 


C. Inclusion of the magnetic field 


For the models discussed above to describe realistic 
vortex lattices magnetic field B must be included in the 
theory. This is achieved by performing the minimal sub¬ 
stitution p — >> p— ^ A in Hamiltonians ([ 5 | and (fir]) where 
A is the vector potential such that B = V x A. Inclusion 
of the magnetic field preserves the discrete symmetries 
listed above but promotes the global U(l) symmetry to 
a gauge symmetry. Specifically, under the transformation 
H(r) e %rZ ^ r ) H{r)e~ lrZx ^ the Hamiltonians remain 
invariant provided that we transform the order parame¬ 
ter phase and the vector potential according to 

0(r) 0(r)-2x(r), (17) 

he 

A(r ) -> A(r) - —Vx(r). 


Here x( r ) is an arbitrary smooth function. It is now 
important to note that while the overlap amplitudes (13) 
and ( fl6| ) are properly invariant under the global U(l) 
symmetry, as written they are not invariant under the 
gauge transformation ( [T7| ). Ref. [51] suggested that in 
the presence of the magnetic field the phase difference 
|(@i — ©j) be replaced by its gauge invariant counterpart 


idij 



•dl, 


(18) 


where 0 j = PkiRj) — 2 Xj(Rj)- The zero mode 

eigenstates are thus given, in this approximation, by Eqs. 
(10) and (15) with Qj replaced by 0j. The overlap in¬ 


tegrals tij in the presence of the magnetic field can be 
computed using Eqs. (13) and (16) with the same re¬ 
placement for Qj. 


An important subtle point here is that Eq. (19) defines 


Xj( r ) only up to an additive constant. The overlap am¬ 
plitudes will have an invariant meaning only if this 
constant is chosen to be the same for all j because then 
it drops out of all differences 0^ — Qj. This condition is 


conveniently implemented by making use of Eq. (18) in 


which the integrand is manifestly gauge invariant, corre¬ 
sponding to a consistent choice of the additive constant. 
Specifically, we conclude that in the presence of magnetic 
field Eq. (16) is replaced by 


t^ = i/iFij sinco’f 


( 22 ) 


with the gauge invariant phase difference defined in Eq. 
(18). The integral itself is path dependent but, as argued 
by Peierls [64], the straight line choice is most physi¬ 
cal because for exponentially localized orbitals the ac¬ 
tual tunneling path is predominantly along the straight 
line where the overlap wavefunction amplitude is max¬ 
imal. The standard Peierls substitution based on this 
reasoning is known to provide an accurate description 
of itinerant electrons moving in ionic lattices subject to 
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magnetic fields. We will demonstrate below, using exten¬ 
sive numerical simulations, that its generalization (18) to 
Majorana fermions in vortex lattices likewise provides a 
description that is both qualitatively and quantitatively 
accurate. 


D. Computation of the phase factors and the Z 2 
gauge structure 


According to the previous subsection in the presence 
of the ma gnet ic field the overlap integrals defined by Eqs. 
(10) and (15) are to be calculated replacing -(0$ — Qj) 


by ujij denned by Eq. ([18]). While the amplitude of ti 3 
depends strongly on various parameters of the model as 
well as on the distance d between the vortices, the Z 2 
gauge factors, which we define as 


Sij = i sgn (sincj^ ) = =bi, 


(23) 


are universal in that they depend only on the vortex lat¬ 
tice geometry. In the following we outline the general 
procedure for the evaluation of these gauge factors and 
we also find them explicitly for some simple lattice ge¬ 
ometries. 

We are mostly interested in physical situations when 
the SC forms a thin quasi-2D layer. In this case the 
effective penetration depth is given by the Pearl length 
A e ff = 2A \/h where \l is the bulk penetration depth and 
h the thickness of the SC film. In most cases we expect 
A l h making A e ff very large. This in turn means that 
the magnetic field can be taken as essentially constant in 
space whenever d < A e ff. 


To calculate the phase factors defined in Eq. (18) 
it is useful to denote the integrand 


n =K ve -l' 1 ) (24) 


and recall that 


V x CL = nz 




R j) $ * 


(25) 


vortex 





FIG. 1: Phase factors and branch cuts in a triangular vor¬ 
tex lattice. Oriented solid lines indicate integration paths 
between the reference points located just to the right of each 
each vortex center. These are used to evaluate the gauge in¬ 
variant phase factors Uij. Dashed lines represent a specific 
choice of the branch cuts discussed in the text. 


of the Cooper pair. Combining Eqs. ( 25|26 ) with the 
Ampere’s law VxB = (47r/c)jf s one obtains the London 
equation for B = zB in the vortex lattice m , 


B-A|V 2 B = $S^(5(r- J R i ), 


(27) 


where \ 2 L = me 2 /47re* 2 n s is the London penetration 
depth. For a periodic lattice the equation can be solved 
by Fourier transforming, 


S(r) = $Si^ T 


iGr 


+ A ?G 2 ’ 


(28) 


where the sum extends over all reciprocal vectors G of 
the vortex lattice. From the knowledge of B one can 
reconstruct n via Eq. ([25]) obtaining 


\ X % iG r 


(29) 


where = he/2e is the SC flux quantum. CL can thus 
be thought of as a vector potential of a fictitious mag¬ 
netic field that consists of ^-function 7 r fluxes associated 
with vortex singularities on top of a neutralizing, almost 
uniform physical magnetic field contributing flux — 7 r per 
vortex. This picture will be useful for determining Uij in 
vortex lattices with high symmetry. One can evaluate CL 
by noting that it is related to the physical supercurrent 
js through 


2/ie* _ , x 

3s=n s — -SI, (26) 

m* 

where n s represents the superfluid density while e* = 2e 
and m* are, respectively, the effective charge and mass 


The gauge invariant phase factors can now be deter¬ 
mined by a straightforward integration of CL(r) in Eq. 
(18) followed by a numerical evaluation of the reciprocal 
lattice vector sums. 

For a regular periodic vortex lattice with high symme¬ 
try, such as the triangular lattice depicted in Fig. [l] it 
is possible to determine the gauge invariant phase fac¬ 
tors uj^ without resorting to detailed calculations. One 
can argue in two stages. First, consider a closed path 
C\ indicated in Fig. [l] It consists of straight line seg¬ 
ments and circular segments. In the following we shall 
consider the latter to have infinitesimal radii. The inte- 
S ral § Ci fl-dl = f(V xfl)- dS can be seen to equal — 37 t/2 
as the path encircles two vortices each contributing flux 
















6 


— 7 r plus an area pierced by magnetic flux 4>o/2 contribut¬ 
ing -b 7 r /2 to the integral. It is also easy to see that the 
circular segments of the path alone contribute the same 
amount of flux — 37 t/ 2. This shows that the total con¬ 
tribution of the straight line segments must be zero. On 
symmetry grounds we furthermore expect each straight 
segment to give the same contribution which must there¬ 
fore be zero. The same conclusion can be reached by 
considering another path, such as C 2 , indicating that this 
result is consistent. We thus arrive at a simple recipe for 
finding ujij\ straight line segments contribute zero while 
the circular segments contribute a phase ±<a/2 where a 
is the angular length of the segment and the sign de¬ 
pends on the sense of rotation with the + sign taken for 
counterclockwise rotation. 

Second, we must attend to the branch cuts. To moti¬ 
vate this consider the path 1 -+ 3 in Fig. [l] According 
to the above recipe we have CJ13 = — 7t/2. However, had 
we avoided the singularity associated with vortex 3 from 
below the result would have been +7t/2. More generally, 
taking the opposite path around the vortex can be seen to 
change ujij —>• ±7r. This ambiguity has to do with the 
fact that as defined in Eq. (24) has an overall factor of 
| in front of the phase gradient which makes sin cjij non¬ 
single valued in the presence of vortices. Importantly, 
this non-single valuedness underlies the Z 2 gauge struc¬ 
ture present in the tight binding model 0 . In order 
to produce a consistent low-energy theory for the MZMs 
we must specify in a globally unique fashion. This 
is achieved by defining branch cuts, emanating one from 
each vortex, along which 17 varies discontinuously. It 
is most convenient to choose the branch cuts such that 
they terminate in a nearby vortex. One such choice of 
the branch cuts is illustrated by dashed lines in Fig. |TJ 
Integration paths between points Ri and Rj chosen so 
as not to cross any branch cuts then furnish a globally 
unique definition of ujij which corresponds to a particular 
choice of the Z 2 gauge. A different choice of the branch 
cuts corresponds to a different Z 2 gauge but leaves all 
physical observables invariant. Factors ±7t/2 indicated 
in Fig. [l]have been obtained according to this prescrip¬ 
tion and can be seen to obey the Grosfeld-Stern rule Eq. 

A similar analysis can be performed for the square 
vortex lattice [5l] and leads to the same conclusion. 


E. Tight binding dispersions for Majorana bands 

We now consider tight binding models for the Majo¬ 
rana zero modes in the triangular and the square vortex 
lattices. The general Hamiltonian is given in Eq. 0 and 
the Z 2 gauge choice is indicated in Fig. [2| Our goal here is 
do derive the corresponding energy dispersions, assuming 
nearest neighbor hopping amplitude t for the triangular 
lattice and both nn and next nn amplitudes t and t' for 
the square lattice. We will then show in the next Sec¬ 
tion that such tight binding models accurately describe 
the MZM dispersions obtained from the full numerical 



b) 



FIG. 2: Vortex lattice geometries: a) square and b) triangu¬ 
lar. Two-vortex unit cell is shaded. The arrows specify the Z 2 
gauge factors for the MZM tight binding models and satisfy 
the Grosfeld-Stern rule Eq. Hopping in the direction of 
the arrow incurs a phase factor of +i while hopping in the 
opposite direction —i. 


solution of the BdG equations describing the p x + ip y 
superconductor and the Fu-Kane model. 


1. Square lattice 

If we denote MZMs associated with the two sublattices 
A and B as qlr and /3r then the Hamiltonian can be 
written as Hq = Hi + H 2 with 

Hi = it olr(Pr - Pn-x-y + Pr~x + Pn-y), (30) 

R 

H 2 = it' ^2 [ a R(~ a R+x + a R+y) + Pr{Pr+x ~ pR+y)\ • 

R 

This can be diagonalized by passing to the Fourier space 



where Q has been inserted for convenience. We note that 
(<+/+ = («—fc-2Q,/?— fc— 2q). If we choose 2 Q = G 
where G is a reciprocal lattice vector we may restrict k 
to one half of the Brillouin zone and regard a^ k and a & as 
regular Dirac fermions defined in the reduced BZ. Differ¬ 
ent choices of Q correspond to different Z 2 gauges. The 
resulting spectra are physically equivalent but may be 
shifted with respect to the center of the BZ. We adjust 
Q as necessary to match the spectra obtained in numeri¬ 
cal simulations discussed below as we do not apriori know 
which gauge is chosen by the numerical diagonalization. 

Choosing Q = (f 5 ~f) allows us to use the 

natural diamond-shaped “antiferromagnetic” BZ for 
this purpose. If we define a two-component spinor 
Tfc = {otk,Pk) T the Hamiltonian takes the form H = 

Efc with 


H k 


fm k h% \ 

V h k ~m k J ’ 


(32) 
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and 


A. Continuum formulation 


hk 


-4 te i ^ kx+ky ^ 2 


sm 


kx T ky 


i sm - 


rrik = 4:t f (cos k x + cos k y ). 


(33) 


The spectrum of excitations 


E k = (34) 

is gapless with a single Dirac point at k = 0 when t' = 0 
and develops a gap A = 8t' otherwise. 


Here we wish to solve the BdG equation (J6| for Hamil¬ 
tonians © and ( fl4| ) defined in the continuum. The 
vortex lattice is encoded in the SC phase field 6{r) as 
described in Eq. 0 and the magnetic field is included 
through the minimal substitution. The key difficulty in 
solving the BdG equation under these conditions lies in 
the fact that although we expect the physical observables 
to exhibit periodicity of the underlying vortex lattice the 
Hamiltonian itself is not periodic. As noted originally 
in Refs. |38l 39 this difficulty can be circumvented by 
performing a singular gauge transformation 


2. Triangular lattice 


H H = UHU~\ 


u = 


/ e -iO A (r ) 

V 0 


o 

e iO B {r ) 


(37) 


For the triangular vortex lattice in the Z 2 gauge indi- where 6U(r) and 0 j g(r) are two functions satisfying 
cated in Fig. [2^b) the MZM Hamiltonian can be written 

as Ha = TL\ + 7 -L 2 with Oa{t) + 0b (r) = 0(r), (38) 


Til = it ^ OLr(—Pr + P R - ai -g 2 - Pr-ch — Pn-g 2 ), 

R 

TL 2 = it + PrPr-clz) . (35) 

R 


Here a± = \/3x and = y are the primitive vectors of 
the sublattice A of the triangular vortex lattice and we 
take the distance between the nn A vortices as our unit of 
length. Fourier transforming according to Eq. (31) with 
the choice Q = (0, 0) leads to the Bloch Hamiltonian of 
the form indicated in Eq. ( |32| with 


and are chosen such that U(r) defined above is single¬ 
valued. In practice this is achieved by partitioning vor¬ 
tices into two sublattices A and B and assigning the con¬ 
tribution from sublattice A to 0a(t) and sublattice B to 
0b{t). The transformed Hamiltonian H is then periodic 
and single valued [38 , 39. and can be analyzed using the 
standard band structure techniques. 


1 . p x + ip y superconductor 

For the p x +ip y SC the transformed Hamiltonian reads 


h k = 4;te^^ kx+ky ^ 
rrik = At sin k y . 


. \/3^r T k y . V^3^x k y 

sm---- + 1 cos-- 


(36) 


The spectrum has the form of Eq. (34) and is fully gapped 


for the triangular lattice. The smallest excitation energy 
4 1 attains at the T point while the maximum 4 y/3t occurs 
along the line between T and M points of the Brillouin 
zone. 


III. NUMERICAL SOLUTIONS OF THE BDG 
PROBLEM WITH VORTICES 

Our objective in this Section is to find numerical so¬ 
lutions of the full BdG equation © in the presence of a 
vortex lattice (and the accompanying magnetic field B ). 
We do this for both the continuum and the lattice for¬ 
mulations of a model p x + ip y superconductor as well as 
the Fu-Kane model. We use these numerical solutions to 
ascertain the validity of the low-energy effective theories 
for the Majorana fermions derived in Sec. II and to relate 
the parameters that enter these theories to the physical 
parameters characterizing the microscopic models. 


H 


= (T(p+<) 2 -» A \ 

V A_ -2 hip- V f) 2 +H)' 


(39) 


where A± = ± iA y and 


A = An 


P+ ~ v f) 


(40) 


Furthermore, quantities 


< = V6»^--A, n = A, B 


(41) 


are related to the physical superfluid velocity v s 

,b 


vf + 

vs and are gauge invariant as well as periodic in real 
space. One can, at least in principle, solve the eigenvalue 
problem defined by H by exploiting the Bloch theorem 
and going to the momentum space. Here, following Read 
and Green [9], we consider a slightly simpler problem 
that follows from sending m 00 in Hamiltonian (|40). 


As argued in Ref. |9] one expects this limit to show the 
same qualitative behavior as the full model: the system 
is in the topological phase with unpaired MZMs in vortex 
cores in the “weak pairing” phase that obtains when /i > 
0 and is in the trivial “strong pairing” phase otherwise. 
















With the above considerations in mind we model vor¬ 
tex core as a small circular region tuned to the trivial 
phase by locally setting (i large and negative. The MZM 
can then be pictured as a chiral edge state at the bound¬ 
ary between the topological bulk and the trivial core re¬ 
gion. This treatment of the vortex core also circumvents 
a difficulty that is known to arise in numerical solutions 
of problems with Dirac Hamiltonians in continuum. As 
first noted in Ref. [39] and later elaborated in Ref. [65] 
eigenstates of the Dirac Hamiltonian tend to diverge as 
~ 1/a/t in the vicinity of the vortex core resulting an 
an overcomplete basis of states. (Note that this behav¬ 
ior is characteristic of the m —>> oo approximation.) To 
treat this problem one must regularize the theory in some 
fashion at short distances. Modeling the core as a triv¬ 
ial strong pairing region represents one possible way to 
regularize by suppressing the wavefunctions at the core 
center. In the next subsection we will discuss the lattice 
formulation of the model which provides another natu¬ 
ral regularization scheme. Importantly, we shall see that 
the low energy properties of the system (i.e. the effective 
theory for the MZMs) are independent of the details of 
the regularization scheme. 

The problem we solve numerically is therefore defined 
by the Hamiltonian 



r X M Y r M 



V A - MO 


(42) 


where A is given in Eq. (40). The superfluid velocities 


that enter the gap function can be determined by a 
procedure analogous to that leading to Eq. (29) above 
(see also Appendix B in Ref. [39 for a more detailed 
description). The chemical potential is taken as 


MO = Mo — M‘ 


'Sr e -(r-Rj) 2 /e 


(43) 


with the second term representing the vortex cores, as 
discussed above. The specific Gaussian form is not im¬ 
portant (any functional form peaked at r = Rj would 
work) but is convenient for the numerics because it has 
a simple Fourier transform. With these preparations we 
can now employ the Bloch theorem, Fourier transform 
the Hamiltonian (42) as described e.g. in Ref. [38], and 
find its energy eigenvalues for each crystal momentum k 
in the first Brillouin zone by a straightforward numerical 
diagonalization. Because the spectrum of the continuum 
Hamiltonian (43) is unbounded we have to impose a high 


energy cutoff A to render the Bloch matrix finite. A must 
be chosen sufficiently large so that the low-energy spec¬ 
trum no longer depends on it. 

Typical results for the square vortex lattice are dis¬ 
played in Fig. [3] In the absence of vortices the spectrum 
shows a gap fi 0 . When vortices are present states ap¬ 
pear inside the gap. These are the expected vortex core 
bound states broadened into bands by intervortex hy¬ 
bridization. The pair of bands closest to zero energy are 
formed of MZMs. We checked that these bands become 


FIG. 3: Band structure for the continuum model of a p x +ip y 
superconductor. In the top panel solid (dotted) lines indicate 
the band structure calculated in the presence (absence) of the 
square vortex lattice. The parameters are chosen as follows, 
= A 0 = 1, p! = 40 and £/a = 0.2. The two bands closest 
to zero energy are the MZM bands. They are enlarged in the 
bottom panel where the dashed line represents the best fit 
to the Majorana tight binding model ( |34| with the hopping 
parameters t — 0.0305 and t' — 0.0076. The inset shows the 
path taken in the first Brillouin zone. 


flat and approach zero energy in the limit of a dilute vor¬ 
tex lattice a £. Their dispersion shows an excellent 
agreement with the tight binding model for MZMs Eq. 
0 with the Z 2 gauge factors given by the Grosfeld-Stern 
rule Eq. (|3|. 


2. Fu-Kane model 


After the singular gauge transformation 
Kane Hamiltonian ( p~4| ) takes the form 


(37) 


the Fu- 


H FK {r) 


fva ■ (p + vf) - p, 

\ A 0 


Ao 

-va ■ (p - vf) 



(44) 

where the vg velocities are given by Eq. (pTb as be¬ 
fore. Once again, the transformed Hamiltonian is pe¬ 
riodic and single valued. Unlike the Hamiltonian for the 
Px + Wy SC which has two distinct phases depending on 
the sign of ji the Fu-Kane Hamiltonian remains in the 
same (topological) phase for all values of // when Aq is 






































9 


non-zero. In order to regularize the wavefunction behav¬ 
ior in the vortex cores we thus introduce a small mod¬ 
ification i^FK -f^FK + 5H m to the Hamiltonian (44), 
making the core magnetic using 


SH r 


•w - C' ra 0 <r) ,Zr)) ■ < 45 > 

As before we take m(r) to be large in the vortex cores 
and zero outside; specifically 


= rn „X^*-(r-R j ) 2 /e 


m(r) = mo > e 


3 


(46) 


and identify £ = v/ttAq with the SC coherence length. 
Magnetic order breaks the time reversal symmetry of 
the TI surface state and is known to gap out the pro¬ 
tected gapless states. The MZMs then can be viewed 
as edge states living on the boundary between the pre¬ 
dominantly magnetic core region and the SC bulk. It is, 
however, important to emphasize that without the long- 
ranged phase structure due to vortices encoded in the v% 
factors the perturbation (45) by itself would not produce 
MZMs as the edge modes would exhibit a large finite size 
gap ~ v H Ao. It is the phase structure that is instru¬ 
mental for the emergence of MZMs whereas m z (r ) serves 
merely to regularize the continuum theory at short dis¬ 
tances. In the next subsection we will see that the same 
MZM structures arise from a theory regularized on the 
lattice where there is no need to include the magnetic 
order. 

With this preparation it is now straightforward to 
numerically diagonalize the Bloch Hamiltonian Hfk{ ft) 
that follows from Eqs. ( 44p6 ) upon Fourier transform¬ 
ing and imposing the high energy cutoff A to render the 
Bloch matrix finite. Typical results of such a calcula¬ 
tion are displayed in Figs. [4] and [5] for the square and 
the triangular vortex lattices, respectively. In both cases 
we observe the initially gapped spectrum (in the absence 
of vortices) modified by the emergence of the low energy 
vortex core states. The bands closest to zero energy arise 
from MZMs. For the square vortex lattice their disper¬ 
sions show near perfect agreement with the Major ana 
tight binding models derived in Sec. II. For the trian¬ 
gular vortex lattice the agreement is also good and can 
be further improved by including longer range hoppings. 
Since the simplest nn tight binding model already cap¬ 
tures all the qualitative features of the MZM band we do 
not pursue this here. 

We conclude this subsection by noting the qualita¬ 
tive and quantitative similarity between the MZM bands 
found in the p x -\-ip y SC and the Fu-Kane model. Indeed 
this is not surprising in view of the expectation that they 
be described by the same minimal tight binding model 
with static Z2 gauge structure described in Sec. II. The 
one distinguishing feature of the Fu-Kane mode - the flat 
MZM bands expected at p = 0 due to the extra chiral 
symmetry - is not apparent in the continuum formula¬ 
tion. This is because the c)i7 m term introduced to regu¬ 
larize the continuum model breaks the chiral symmetry 




FIG. 4: Band structure for the continuum Fu-Kane model, 
square vortex lattice. In the top panel solid (dashed) lines in¬ 
dicate the band structure calculated in the presence (absence) 
of the vortex lattice. The parameters are chosen as follows: 
v = 0.14, mo = 7.0, p = 0 and £/a = 0.14. The high energy 
cutoff A = 8 and all the quantities are in units of Ao = 1. 
As in Fig. [4] the two bands closest to zero energy are the 
MZM bands. They are enlarged in the bottom panel where 
the dashed line represents the best fit to the Majorana tight 
binding model ( p34| with the hopping parameters t — 0.0153 
and t! = 0.0032. 


(even at p = 0). We were unable to find a symmetry pre¬ 
serving regulator that would work for this purpose in the 
continuum model. We shall see however that the lattice 
model considered next preserves the chiral symmetry and 
indeed exhibits the expected flat bands at zero energy. 


B. Lattice formulation 

Although technically somewhat more complicated the 
lattice formulation of the problem has a distinct advan¬ 
tage of providing a natural short distance cutoff for the 
electron wavefunctions in vortex cores. Artificial regu¬ 
lators that were necessary in the continuum theory are 
thus not needed. The lattice formulation of the p x + ip y 
superconductor with vortices has been discussed in Ref. 
[39] although Majorana bands have not been studied in 
detail. Here we briefly review the lattice construction 
and examine the Majorana bands. This we follow by a 
similar discussion for the Fu-Kane model whose lattice 
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FIG. 5: Band structure for the continuum Fu-Kane model, 
triangular vortex lattice. Solid (dotted) lines indicate the 
band structure calculated in the presence (absence) of the 
vortex lattice. The parameters are chosen as follows: v = 
0.33, m = 1.0, fi — 0 and £/a = 0.17. The high energy cutoff 
A = 20 and all the quantities are in units of Ao = 1. As in Fig. 
[4] the two bands closest to zero energy are the MZM bands. 
The dashed lin e rep resents the best fit to the Majorana tight 
binding model (35) with the hopping parameter t — 0.0425. 


formulation has not been previously discussed. 


1 . p x + ip y superconductor 


The problem is defined by the Hamiltonian (4|5 ) with 
the fermion creation operators now residing on sites of a 
lattice which we take to be square and of unit spacing. 
The kinetic and pairing operators are given by [39 


h 

= ~rJ2 e~ i(e/hc) K +S A ( r > dl s s - e F , 

§ 

(47) 

A 

= A 0 Y,e ie{r)/2 fise i9ir)/2 . 

(48) 


5 


From now on, we will set the hopping amplitude r to 
unity and measure all energies in units of r. Also, s$ 
denotes the shift operator sgu(r ) = u(r + 5), where 5 
represents a nn vector. For the p x + ip y superconductor 
the operator i)$ is defined as 


Vs 


Tiss if 5 = d zx 7 
±ss if S = ±y. 


(49) 


After the singular gauge transformation (37) we obtain 
the Hamiltonian 


(-Ee iV ° ir) s s -e F 

H = s 

\ <5 

with the phase factors defined as 

v£( (r) = C +S ( W/i “ ¥c A ) ' dl v = A ’ B ( 51 ) 


Vs 


Ao£e^ } 

'£e- i v?Ws6+e F 

s 


(50) 




FIG. 6: Band structure for the lattice Px+Wy superconductor 
with a square vortex lattice, a) Full band structure in a 10 X10 
magnetic unit cell with Ao = 0.50 and £f — —2.2. b) Detail 
of the Majorana band (solid squares) and the best fit to the 
tight-binding Majorana dispersion Eq. (34) with t — 4.3027 x 
10 -3 and t! — 1.9375 x 10 -4 (solid line). 


and As(r) = |[V^(r) — V$(r)]. The phase factors are 
easily evaluated using the method discussed in Sec. II.D 
(see also Appendix B in Ref. [39] for details). 

The Hamiltonian ( |50| ) now has the periodicity of the 
vortex lattice (with two vortices per unit cell) and can 
be diagonalized in momentum space using standard band 
structure techniques. Fig.|6ja) shows the band structure 
obtained for the square vortex lattice. It exhibits the 
expected Majorana band close to zero energy as well as 
nearly flat Landau level bands at energies high compared 
to the SC gap amplitude Ao, in complete agreement with 
results of Ref. [39] . Panel (b) of the figure focuses on the 
Majorana band which is, once again, very well described 
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weak coupling limit A 0 <C Sp with Sp referenced to the 
bottom of the band. Here R denotes the distance between 
the vortices, kp is the Fermi momentum and £ = vf/ttAq 
is the BCS coherence length. The above expression ac¬ 
curately captures the period and the phase of oscillations 
in both t and t' but does not describe the amplitude par¬ 
ticularly well. We tried other, more complicated expres¬ 
sions derived in Ref. [48} , but they do not significantly 
improve the agreement. We instead find that the data is 
well described by a simple phenomenological expression 


t ~ 4(1 + bRk) | cos (kR + 7r/4)|. (53) 

Here k = \Jk\ — ( Aq/vp ) 2 while A and b are dimen¬ 
sionless constants, and a similar expression for t' with 
parameters A' and b'. Because the MZM wavefunctions 
decay exponentially with the characteristic lengthscale 
£ the amplitude implied by the expression (52) makes 


good intuitive sense. Our results suggest that the inter¬ 
play between MZM wavefunctions in the vortex lattice is 
possibly quite intricate and cannot be fully captured by 
the perturbative treatment of two distant vortices carried 
out in Refs. 03 HE]. 


FIG. 7: The Majorana overlap amplitudes t and t' as a func¬ 
tion of £f extracted from the lattice model p x + ip y super¬ 
conductor (solid symbols). The magnetic unit cell is 50 x 50 
and Ao = 0.50. The thin (blue) line represents the analytical 
result of Ref. [48] while the thick (red) line corresponds to 
the simple phenomenological expression ( [53] ) discussed in the 
text. The fit parameters are A = 2.43 x 10 -4 , b = 0.0163 and 
A' = 5.91 x FT 5 , b' = 0.0066. 


by the tight-binding dispersion (34). 


From these results we may easily extract the depen¬ 
dence of the tunneling amplitudes t and t' that enter the 
effective tight-binding Majorana model on various micro¬ 
scopic parameters of the underlying BdG theory as well 
as the vortex lattice geometry. We do this by fitting the 
numerically calculated MZM bands to the tight binding 
dispersion ( [34] ) . As an example Fig. [7] displays the de¬ 
pendence of ft, t') on the Fermi energy £p with Aq, t and 
vortex spacing held fixed. We note that in view of Eq. 
(34) the Majorana band structure is not sensitive to the 


sign of the amplitudes t and t'. We may plausibly surmise 
that the nodes apparent in Fig. [7] represent sign changes 
in the amplitudes which fixes them up to an overall sign. 
The overall sign could be potentially determined from 
the structure of the corresponding wavefunctions but we 
do not pursue this issue here since we do not believe the 
sign is an easily measurable quantity. 

We can compare these with the analytical expressions 
derived in Refs. mm- Specifically, we plot 


2 \cos (k F R+ J)| 

— vm — expl 


(-f) <»> 

which corresponds to Eq. (31) of Ref. [48], valid in the 


2. Fu-Kane model 

The implementation of the Fu-Kane model on the lat¬ 
tice is more involved owing to the Nielsen-Ninomyia the¬ 
orem [66] . which states that it is impossible, as a matter 
of principle, to construct a T-invariant 2D lattice Hamil¬ 
tonian with an odd number of Dirac fermions in the low- 
energy spectrum. It is therefore impossible to write a 2D 
lattice model that would faithfully describe a single sur¬ 
face of a TI. Studying the full 3D problem (including the 
STI bulk) would provide the desired outcome but this 
would be computationally very costly. We also note that 
magnetic field of several Tesla, sufficient to produce the 
vortex lattice, has negligible effect on the gapped bulk of 
the STI. This is because the relevant cyclotron frequency 
as well as the Zeemann energy are much smaller than the 
bandgap (which is ^ 300 meV in B^Ses family of mate¬ 
rials). There is, therefore, nothing interesting to learn by 
performing a full 3D calculation. The problem of doped 
STI, where the bulk itself can become superconducting, 
has been studied with some interesting results [67,1.68J. 

To circumvent the above problem we employ the idea 
introduced in Ref. [69] and construct a lattice model de¬ 
scribing instead a pair of parallel TI surfaces, such as 
those terminating a slab. Because a pair of TI surfaces 
has in general an even number of Dirac fermions the the¬ 
orem 66 no longer presents an obstruction. Ref. [691 
showed how to construct a lattice model of this type with 
low-energy degrees of freedom on two surfaces that are 
largely decoupled. 

The normal state Hamiltonian (“model II” in Ref. 69]5 
can be written in the momentum space as 

r _ ( g_k M k \ 
k ~ \M k -g k ) 


( 54 ) 
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with = 2X(a y sin k x — a x smk y ) and = 2r(2 — 
cos k x — cosky). Its diagonal blocks describe the gapless 
surface states in the two surfaces of a TI and the coupling 
Mk is designed to gap out all the Dirac nodes except 
those at the origin, k = (0,0). In the following we set 
A = 1 and measure all energies in units of A. To study the 
vortex lattice we imagine inducing superconductivity in 
one of the surfaces by proximity effect. This is described 
by passing to the BdG formulation using Eq. (14) with 


A = 


A 0 
0 0 


(55) 


To avoid complications that would arise from the other 
surface being ungapped we imagine that its surface state 
has been gapped by a T-breaking perturbation and re¬ 
place gk Qk + rna z in the lower diagonal element of 
Eq. (54). Since the physical time-reversal is already bro¬ 


ken by the applied magnetic field we do not expect this 
additional T-breaking to have a significant effect on the 
system, other than removing the unwanted gapless exci¬ 
tations from the second surface. We also emphasize that 
this is a purely technical device and we do not require 
such magnetization to be implemented in the experimen¬ 
tal realization. 

The full momentum-space Hamiltonian we consider 
thus has the following form 


K K 


(gk — £f 
M k 
A* 

0 


M k 

~g k ~ ma z 
0 
0 


A 

0 

~9k + £f 

-Mk 


\ 


o 
o_ 

-Mk 

g k ~ ma z J 
(56) 

In order to implement the vortex lattice we now pass to 
the real space and perform the minimal substitution to 
include the magnetic field. The upper diagonal block of 
H fk thus becomes 


/ Sf 


s 

5 

4 t-t^2§s 

s 

0 


4 t-t'Ess 0 \ 

—sf 0 4 t — t^ss 

s 

0 —m 

d 


V 


4 t-t^ss ~iJ2vs 

s d 


/ 


and a similar expression for the lower diagonal block. 
Vortices are included by replacing A —>* Ae l6 ^ and the 
magnetic field enters via the Peierls substitution 


ss 


0 -i(e/hc ) f^ +S A(r)-dl 


Ss • 


(57) 
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FIG. 8: Band structure for the lattice version of the Fu-Kane 


model with a square vortex lattice, a) Full band structure in 
a 30 x 30 magnetic unit cell with r = 0.5, Ao = 0.4, m = 0.5 
and £f — 0.25. b) Detail of the Majorana band (solid squares) 


and the best fit to the tight-binding Majorana dispersion Eq. 
(34) with t = 2.36 x 10 -3 and t' = 3.45 x 10 -4 (solid line). 


The band structure for a square vortex lattice and a 
generic Fermi energy sp is displayed in Fig. [8j The Ma¬ 
jorana band shows a weak dispersion and is, once again, 
well described by the effective Majorana tight binding 
model discussed in Sec. II.E. One can extract the tunnel¬ 
ing amplitudes t and these are plotted in Fig. [9] We 
find that they are reasonably well described by a simple 
heuristic formula 


We note that, importantly, the Peierls phase factors must 
now be also attached to the shift operators that enter 
the definition of i)s because in the Fu-Kane model these 
appear in the kinetic energy of the system and thus rep¬ 
resent single-electron hopping processes. 

As before, singular gauge transformation (37) renders 
the Hamiltonian periodic and we can solve it in mo¬ 
mentum space using standard band structure techniques. 


t Ae bkirR |sin (kpR) \ , (58) 

and a similar expression for t' with parameters A' and b '. 

When sp is tuned to the neutrality point we observe 
that both t and t' vanish which results in a completely 
flat Majorana band, as expected in the presence of the 
extra chiral symmetry discussed in Sec. II.B. While the 
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r x M Y r m 


FIG. 10: Band structure for the lattice version of the Fu-Kane 
model with a deformed square vortex lattice, as explained in 
the text. A 30 x 30 magnetic unit cell is used with vortices 
located at (5,5) and (-5,-5) basis vectors. Also r = 0.5, Ao = 
0.4, rri — 0.5 and ef — 0.0. 


IV. CONCLUSIONS 


FIG. 9: The Majorana overlap amplitudes t and t' as a func¬ 
tion of ef extracted from the Fu-Kane model formulated on 
the lattice (solid symbols). The magnetic unit cell is 30 x 30, 
r =s 0.5, Ao a=s 0.4 and m — 0.5. The solid (red) line cor¬ 
responds to the simple phenomenological expression (57) dis¬ 
cussed in the text. The fit parameters are A = 7.85 x 10 —3 , 
b = 0.224 and A' = 5.87 x 10 -4 , b' = 0.170. 


tunneling amplitudes also vanish for certain nonzero val¬ 
ues of £p, these are accidental zeros. Importantly, £p = 0 
is the only value for which t and t' (and presumably all 
other amplitudes) vanish simultaneously. To test that 
the flat band is indeed protected by the chiral symmetry 
(and not by some symmetry of the square vortex lattice) 
we performed simulations for a “deformed” square lat¬ 
tice. It is defined as follows: we keep the unit cell the 
same but within the unit cell we gradually move the B 
vortex closer to the A vortex along the diagonal line that 
connects them. When £f — 0 the Majorana band re¬ 
mains completely flat for each A-B vortex separation, all 
the way to the point when the two vortices merge and 
form a square lattice of doubly quantized vortices. This 
is illustrated in Fig. [lO] 

i 


Majorana zero modes bound to vortices in topological 
superconductors form bands in the presence of a vortex 
lattice. We have demonstrated that such bands are well 
described by simple tight binding models Eq. 0 describ¬ 
ing short ranged tunneling events between the adjacent 
sites on the lattice. An interesting feature of these mod¬ 
els is the underlying non-trivial Z 2 gauge structure that 
is mandated by the canonical anticommutation relations 
for the self-adjoint MZM operators §. We found that 
when the magnetic field necessary for the vortex lattice 
formation is properly included the Z 2 gauge factors obey 
the Grosfeld-Stern rule 0 , previously derived in the con¬ 
text of MZMs in the Moore-Read fractional quantum Hall 
state as well as vortex lattices mam] in the Kitaev spin 
model on the honeycomb lattice [45]. The hopping am¬ 
plitudes are found to retain the full periodicity of the 
vortex lattice and do not show any anomalies suggested 
by previous works that neglect the applied magnetic field 
[33H35]. 

For periodic vortex lattices, such as the square and 
the triangular lattice, the resulting low-energy theory is 
typically gapped and topologically nontrivial. The latter 
property follows from the non-zero gauge flux implied by 
the Grosfeld-Stern phase factors. These, in turn, orig¬ 
inate from the structure of the individual MZM wave- 
functions and their overlap integrals. An intuitive under¬ 
standing of this structure can be obtained from the fol- 
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lowing simple argument. Because of the self-adjoint prop¬ 
erty ([2]) of the MZM operators and their fermionic anti¬ 
commutation relations the hopping amplitude between 
two sites is necessarily imaginary. This corresponds to 
the Z 2 phase ±i. Now the simplest closed path on the 
lattice involves 3 distinct sites. The total Z 2 phase ac¬ 
cumulated along such path is ±i which corresponds to 
a non-zero enclosed Z 2 flux. It follows that, generically, 
Majorana fermions defined on a lattice move in the back¬ 
ground of a non-vanishing Z 2 gauge flux. In analogy 
with the Haldane model 12a one then expects the system 
to exhibit a non-zero Chern number and, in the geome¬ 
try with open boundaries, protected gapless edge modes. 
This indeed has been noted in previous theoretical stud¬ 
ies [34] [35]. 

In the Fu-Kane model an interesting situation arises 
near the so called neutrality point where an extra chiral 
symmetry exists. The latter mandates that all hopping 
amplitudes exactly vanish resulting in the MZM band 
that is completely flat. This expectation is indeed borne 
out by our analytical as well as numerical calculations. 
Such completely flat bands are then highly susceptible to 
the effects of interactions and disorder. Some of the inter¬ 
action and disorder effects have been explored in recent 
works [5T1 [56] and found various interesting interacting 
phases of Majorana zero modes in one and two dimen¬ 
sions. 

In a homogeneous superconductor the vortex lattice 
is expected to be perfectly periodic GE3GE] and our re¬ 
sults then directly apply. Many clean superconductors 
indeed exhibit such perfectly periodic vortex lattices. In 
a disordered superconductor, however, vortex lattice it¬ 
self may become disordered. Our method for calculating 
the full electronic structure relies on translational invari¬ 
ance and cannot be directly applied to such disordered 
vortex lattices. However, our results indicate that the 
effective Majorana tight-binding model Eq. 0 provides 
a good description of the low-energy physics. One thus 
expects that Majorana degrees of freedom in a disordered 
vortex lattice will be well described by the same tight- 
binding model in which the overlap integrals acquire a 
random component. Randomness in this model has been 


extensively studied EMU and we expect these results to 
directly transfer to the present problem of vortex lattices 
with randomness. In addition, in a 2D system individ¬ 
ual vortices can undergo thermal or quantum fluctua¬ 
tions around their equilibrium positions. The fate of the 
energy bands in this situation is an interesting problem 
which we leave for future study. 

Recently, individual vortices have been experimentally 
observed by scanning tunneling microscopy (STM) in 
2D heterostructures combining a topological insulator 
Bi 2 Te3 and a conventional superconductor NbSe 2 {28] . 
Evidence for possible MZMs bound in the cores of such 
vortices has also been reported [29] . Although these ex¬ 
periments were not performed in the parameter regime 
where the MZM band formation could be directly ob¬ 
served these developments suggest that the results ob¬ 
tained in the present work can be experimentally tested 
in the near future. 

After this work was submitted for publication we be¬ 
came aware of a preprint m that reports results on the 
band structure of p x + ip y superconductor with vortices 
in agreement with our results. The preprint also studies 
in detail the topological phases of the Majorana bands 
and finds results that support conjectures presented in 
this section, in addition to many other interesting new 
results. 
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